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We introduce a new class of two(multi)-matrix models of positive Hermitean matrices coupled in a chain; 
the coupling is related to the Cauchy kernel and differs from the exponential coupling more commonly 
used in similar models. The correlation functions are expressed entirely in terms of certain biorthogonal 
polynomials and solutions of appropriate Riemann-Hilbert problems, thus paving the way to a steepest 
descent analysis and universality results. The interpretation of the formal expansion of the partition 
function in terms of multicolored ribbon-graphs is provided and a connection to the O(l) model. A 
steepest descent analysis of the partition function reveals that the model is related to a trigonal curve 



(three-sheeted covering of the plane) much in the same way as the Hermitean matrix model is related to 
a hyperelliptic curve. 
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1 Introduction 

In the last two decades or so, the theory of matrix models have been an incredibly fertile ground for a 
fruitful interaction between theoretical physics, statistics, analysis, number theory and dynamical systems 
(see the classical [35] and references therein). 

The interplay with analysis has been particularly beneficial for the Hermitean matrix model [14], 
largely due to the realization that the matrix model could be "solved" in terms of orthogonal polynomials. 
This allowed to translate questions about the spectrum of large random matrices into questions about 
the asymptotics of orthogonal polynomials, connecting the former with a very well developed area of 
analysis. The availability of Riemann-Hilbert methods [28, 13] was the crucial ingredient in addressing 
questions of universality in the bulk and at the edge of the spectrum. 

It should be mentioned that this fortunate symbiosis relies on the following features: 

• the possibility of rewriting the matrix integral in terms of eigenvalues and the correlation functions 
in terms of suitable orthogonal polynomials; 

• the Riemann-Hilbert characterization of orthogonal polynomials [28] ; 

• the Christoffcl Darboux formula, allowing to rewrite the kernel of the correlation functions in terms 
of only two orthogonal polynomials of degree TV; 

• the nonlinear steepest descent method [13] applied to the RH problem for orthogonal polynomials. 



There are by now several matrix models based on various ensembles of matrices. For some (Symplectic, 
Orthogonal), a connection with (skew-) orthogonal polynomials can be established [35, 37] and some of the 
steps in the above list have been performed, but typically only for a certain restricted class of "potentials" . 
For others (like 0(n) models, rectangular models) the methods used rely on symmetries of the integral 
and dynamical system approaches (the "loop equations") [18, 16, 17, 20, 1] 

There are also the so-called multi-matrix-models which involve ensembles of several matrices (typically 
Hermitean) and one of the most studied among them is the Itzykson-Zuber-Harish-Chandra (IZHC) chain 
of matrices. In its simplest form, the two-matrix model, it amounts to the study of the spectral properties 
of a pair of Hermitean matrices of size N x N with a measure 

d l x(M 1 ,M 2 ) = dM 1 dM a e- NTt<y ^ Ml)+v ^ M ^- MlM ^ (1-1) 

where the interaction term is e JVTrM i M 2 p xhe first three items of the bulleted list above can be imple- 
mented at least for potentials whose derivative is a rational function [19, 6, 4, 5] while the last item is still 
not fully under control. One of the main reasons for this difficulty is that the size of the RHP depends 
on the potentials: for example if Vj are polynomials of degree dj then there are two relevant RHPs of 
the same size. Clearly the question arises as to whether a Ricmann Hilbert method can be utilized for a 
more general class of potentials, for example real- analytic, like in the case of the single-matrix model. 

The present paper is a part of a larger project initiated in [8], which puts forward a new multimatrix 
model (at the moment primarily two matrices) that can be completely solved along the lines described 
above. 

While the new model is very close in spirit to the IZHC model, the different interaction 

^^ det^ + M,)" M 

links the model to a new class of biorthogonal polynomials which were studied in [8] and termed "Cauchy 
biorthogonal polynomials" . These polynomials have several of the desirable features of classical orthogonal 
polynomials. 

We will fill in more details in Sect. 3 but here we just indicate that the model we want to study is 
defined on the space of pairs of positive Hermitean matrices of size N equipped with a measure 

for arbitrary measures a(x)dx,(3(y)dy on E + (to be understood in the formula above as conjugation- 
invariant measures on positive definite Hermitean matrices or, equivalently, measures on the spectra of 
Mi's). This immediately puts the problem at the same level of generality as the classical case. Cor- 
responding to the above positive measure is the normalizing factor, customarily called the partition 
function 

Z N := f I 'd Ml dM 2 «j, (1-4) 

J J 2 det(Af! + M 2 ) w l ' 



whose dependence on the measures a,/3 carries all relevant information about the model. 
We will show that this model is related to Cauchy biorthogonal polynomials 

a(x)dx (3(y)dy 

Pn{xjq m (yj , = <W (1-5) 

x + y 

defined and studied in [8] in relation with the spectral theory of the cubic string and the Degasperis- 
Procesi wave equation (see also [2, 3, 34, 33]). Contrary to the biorthogonal polynomials of [19, 6, 7] the 
properties of these polynomials do not depend on the measures a, (3 and the main highlights are 

1. they solve a four-term recurrence relation; 

2. their zeroes are positive and simple; 

3. their zeroes have the interlacing property; 

4. they possess Christoffel-Darboux identities relevant to matrix models; 

5. they can be characterized by a (pair of) Riemann-Hilbert problem(s) of size 3x3; 

6. the steepest descent method is fully applicable. 

Our recent paper [8] contains a detailed discussion of all these points except the last one which will be 
addressed in a forthcoming publication. 

The paper is organized as follows; in Sect. 2 we review our previous results [8] on the Cauchy 
biorthogonal polynomials and the relevant formulae and features needed in the following. In particular 
the Christoffcll-Darboux identities (Sec. 2.1), the Riemann-Hilbert characterization in terms of 3 x 3 
piecewise analytic matrices (Sect. 2.2). In Prop. 2.7 we introduce the 3x3 matrix kernel in terms of 
the solution of the Riemann-Hilbert problems; this will be used in later section to describe the spectral 
statistics of the two matrices. In Sect. 3 we introduce in detail the two-matrix model that we have outlined 
above and show how to reduce the study of its spectral statistics to Cauchy biorthogonal polynomials 
using a formula appeared in [25]. We also indicate (without details) how to deal with a similar model 
where R matrices are linked in a chain (Sect. 3.2). 

In Sect. 4 we show how a formal treatment of the partition function of the model for large sizes of 
the matrices can be used to extract combinatorial information for certain bi-colored ribbon graphs (as 
it has been done for the Hcrmitean one-matrix model in [11, 32] and for the multi-matrix models with 
exponential coupling in [16]). 

In Sect. 5 we show that a saddle-point treatment of the partition function (see for example [15] for the 
one-matrix case) leads to a three-sheeted covering of the spectral plane (a trigonal curve) ; this (pseudo) 
algebraic curve plays the same role as the hyperelliptic curve in the Hermitean matrix model. 

The Appendices are devoted to the relation between our proposed model and other matrix models 
with rectangular matrices (and possibly with Grassmann entries, App. A) and a connection (App. B 
with the O(l) model of self-avoiding loops ([18] and references therein). 



Remark 1.1 For a particular case of measures a{x) = x a e~ x and f3(y) = y b e~ y the corresponding 
biorthogonal polynomials appeared (in a somewhat disguised form) in the work [10]. As observed therein 
they are related to the classical Jacobi orthogonal polynomials for the weight x a+b dx on [0, 1]. We thank 
A. Borodin for pointing out this connection. 

2 Cauchy biorthogonal polynomials 

Let K(x,y) — -^r- be the Cauchy kernel on M. + x R + . It is known that it is totally positive in the 
sense of the classical definition [22, 31]: 

Definition 2.1 A totally positive kernel K(x,y) on 3 x $ c M. x M. is a function such that for all 
n € N and ordered n-tuples x\ < X2 < • • • < x n , y\ < y% < ■ ■ ■ < y n we have the strict inequality 

det[K( Xi , yj )}>0. (2-1) 

It was shown in [8] that for any totally positive kernel K(x, y) (hence also for the Cauchy kernel) and for 
any pair of measures a(x)dx,/3(y)dy supported in 3,3 C R+ respectively, the matrix of bimoments 



hj = ^ J^ x % y>K{x, y)a(x)dx(3(y)dy (2-2) 

is a totally positive matrix, that is, every square submatrix has a positive determinant. This guarantees 
the existence of biorthogonal polynomials {pj(x), Qj(y)} of exact degree j such that 



Pj(x)qk(y)K(x,y)a(x)dx(3(y)dy = 5 ]k . (2-3) 

The polynomials pj(x),qj(y) are defined uniquely up to a C x action pj t— > XjPj(x),qj i— ► j^lj & n d the 
ambiguity can be disposed of by requiring that their leading coefficients are the same and positive. With 
this understanding we have p n (x) = —x n + . . . and q n — —y n + . . . , with the positive constant c n given 
by [8] 

c » = xl^rT 1 ' Dn := det [ J ^]o<,-fc<n-l > ^ ■= [ [ ^y k K{x 1 y)a{x)dxP{y)dy . (2-4) 

V J-'n - J ' - Jr + Jr + 

The determinantal expressions for the BOPs in terms of bimoments can be obtained using Cramer's 
rule [8]. 

Proposition 2.1 (Thm. 4.5 in [8]) For any totally positive kernel K(x,y) the zeroes of the biorthogo- 
nal polynomials pj 's (qj 's) are simple, real and contained in the convex hull of the support of the measure 
a (P respectively). 

In the case K(x,y) = -^— we named the corresponding polynomials {pj(x), qj(y)} Cauchy BOPs 
and proved, in addition, 



Proposition 2.2 (Thm. 5.2 in [8]) The roots of adjacent polynomials in the sequences {pj(x)}, {qj(y)} 
are interlaced. 

Cauchy BOPs enjoy more structure than the BOPs associated to a generic kernel: they solve a four 
term recurrence relation of the form 

x(x n -iPn{x) - it n p n -i{x)) = c4 -1) pn+i(a;) + a^p n (x) + a£>p n -i(x) + a^'p n -2(x) 
y(vn-iq n (y) - Vnq n -i(y)) = ^ _1) 9n+i(y) + b^q n (y) + b^qn-iiy) + b^ q n _ 2 (y) (2-5) 

where 

7T„ := / p n (x)a(x)dx , i] n := q n (y)[3(y)dy . (2-6) 

As proved in [8] n n ,ri n are strictly positive. The other coefficients appearing in the recurrence relations 
are described in loc. cit. 

We also need to introduce certain auxiliary polynomials, q n ,p n whose defining properties are 

1. degcf n = n + l, degp„ = n; 

2 I q n &(3 = or q n (y) = ^^ " ^ 



3. // Pn(x)q m (y) ; = S mn Or -^-(Pn-l ~ Pn) = Pn- 

J J x + y *>n 

In addition p n , q n admit the determinantal representations: 
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VnVn+l\/ 'D n D n+2 
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(2-7) 



(2-8) 



•• Pn 

After lengthy manipulations one obtains several Christoffel-Darboux-like identities (CDIs) for Cauchy 
BOPs which play a crucial role in what follows and hence will be carefully described. 

2.1 Christoffel— Darboux Identities 

Using recurrence coefficients featured in eqs. (2-5) we define the following 3x3 matrices 



A n (x) = 









b {2) 1 

°n+l 


, »n(y) = 








(2) " 
a n+l 


°n-l 



° n + rin+i 

-b { - 1] 






-a (_1) 




„(0) . 3/ 






fl " ' TTn + l 

-o (_1) 



(2-9) 



In addition we define the following integral transforms 



W(*\ ■= 



Pi L, (z) 



and the following vectors 



q^(z) = 



w — C 
Pn(Qad£ 






Mi 



(11/ v /* 9«(C)/?dC ( 2 )/ s f 9n'(-X) 



10 + X 



-ct(x)da; 



1+ z + V 



(2-10) 
(2-11) 



p£ m) (*) = 



P^i(z) 
p£°(s) 



(2-12) 



where /^ = 0,1,2 and g„ = g n , p„ = p n . For a(x),{3(y) we define a* (a;) = a(—x) and /3*(y) = (3(—y); 
next we define the Weyl functions (or Markov-functions) as 

1 



W fi (z) 



z-y 



0{y)dy = -Wp*(-z), 



W a ,{z) 



-a{x)dx — —W a (—z), 



W a , p {z) = - jj ^ a(x)/3(y)dxdy W 0a ,(z) = jj * a(x)P(y)dxdy (2-13) 

Proposition 2.3 (Thm 7.3 in [8]) The following Christoffel-Darboux-like identities hold 

n-i 
(z + w)J2 q^\w)pf\z) = qW^ • A(-ta) • K M (z) - ¥(w, *)„,* (2-14) 

j=o 

"0 1 

¥(w,z)= 1 Wi3*(z)+Wp(w) (2-15) 

1 W„(«) + W a .(«0 W„.(«7)W / J.(«)+Wa^(w) + W / J. a (2). 

where the auxiliary vectors marked with a hat are characterized by a Riemann-Hilbert problem described 
below. 



Corollary 2.1 (Thm 7.4 in [8]) Evaluating (2-14) on the "antidiagonal" z 
duality 



1 

1 

1 



-w gives the perfect 



(2-16) 



Remark 2.1 The proposition above defines, in fact, 9 identities, but we will only need the 4 identies 
corresponding to /U, v = 0, 1 (i.e. the principal submatrix of size 2 x 2). 



2.2 Riemann— Hilbert characterization of the integrable kernels 

The sums appearing on the left hand side in Prop. 2.3 are all examples a general framework of "integrable 
kernels" that were studied in great generality in [30, 27, 29, 26, 24]. 

Proposition 2.4 ( Prop. 8.1 in [8]) Consider the Riemann-Hilbert problem (RHP) of finding a ma- 
trix T(w) such that 

1. T(w) is analytic on C\ (supp(P) U supp{a*)) 

2. T(w) satisfies the jump conditions 



r(«>) + = r(u>). 
r(«>) + = r(t«). 



1 -2nif3 
1 

1 

1 

1 -2-Kia* 

1 



w G supp{(3) C R_| 
w G supp(a*) C - 



3. its asymptotic behavior at w — oo Q(w) ^ is 

T(w) = (1 + Oiur 1 )) 



w n 
w- 1 

w- n+1 



Then such a T(w) exists and is unique. Moreover T(w) can equivalently be written as: 



or also 



r(«0 
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y(w) 



J(l) 



(2)/ 



(2-17) 



(2-18) 



(2-19) 



q i 0) M , q W (w) , qy } (w)] (2-20) 



where the normalization constants r/ n ,c n have been introduced in (2-4, 2-6). 

Remark 2.2 The quantities referring to the letter q on the right hand side in Prop. 2.3 can be extracted 
from the RHP involving T. The next proposition will achieve the same goal for the remaining quantities. 

Proposition 2.5 (Prop. 8.2 in [8]) Consider the Riemann-Hilbert problem (RHP) of finding a matrix 
r(z) such that 



1. r(z) is analytic on C \ {supp(a) U supp{f3*)) 



2. T(z) 


satisfies the jump conditions 












1 -2ma(z) " 






r(z) + = r(z)_ 


1 
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, z € supp(a) C K + 




f(z) + = f(*)_ 
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1 -2ni(3* 


, z e supp(/3*) CR_, 
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(2-21) 



3. its asymptotic behavior at z — oo ^s(z) ^ is 



T(z) 



O 



z n 
1 
ir 



(2-22) 



TTien smc/i a T(z) exists and is unique. Moreover T(z) can equivalently be written as: 

T(z) 



c n 

0-10 

o o ^ 



P0,n Pl,n P2,n 

P0,n-1 Pl,n-1 P2,n-1 

P0,n-1 Pl,n-1 P2,n-1 



(2-23) 



or also 



T{z) 







(-1)' 





-1 










r(z) 



1 -1 o 
1 
-1 1 



AO) zfl) ^(2) 

P„ (z),P„ {z),P n {z) 



(2-24) 



c„_i?7 n _i 

where the normalization constants rj ni c n have been introduced in (2-4, 2-6). 

Remark 2.3 These polynomials p n {z) :— p\ '(z) and the auxiliary functions ph (z),pn (z) were intro- 
duced in [8] independently from a RHP formulation, but for all practical purposes the formulation above 
is sufficiently explicit. 

Finally the RHPs above allow us to reconstruct the ratio of two consecutive principal minors of the 
bimomcnt matrix (this will become relevant when discussing the partition function of the matrix model) 

Proposition 2.6 (Corollary 8.1 in [8]) If D n is a leading principal n x n minor of the bimoment 
matrix I , then 



D„ 



= (-1)" lim w< 



% -iT%,i{w) 



<-i > o 



(2-25) 



D n _i "' «J-»"5o r 2 ,i(w) 

The CDIs in Prop. 2.3 can be written in a very simple form in terms of the solutions of the Ricmann- 
Hilbert problems. To show this let us momentarily denote by Y(w) and Y(z) the matrices with columns 
given by the q„ (w) and p„ (z) respectively. Then (2.3) can be reformulated as 



(z + w)M n (z, w) :— {z + w) 



n-l 

J2#\w)pf\z 

3=0 



+ ¥(w, z) = Y\w)K{-w) ■ Y{z) 



(2-26) 



Now the perfect duality of Cor. 2.1 implies that 

Y(w)A(-w) = W(-w)- 1 . (2-27) 

Note that the solutions of the RHPs T(w),T(z) differ from Y(w),Y(z) only by some constant invertible 
left multipliers, and hence 

Y(w) ■ A(-to) • Y(z) = J • Y(-w)- 1 ■ Y{z) = J • f (-w)" 1 • T(z) . (2-28) 

We collect this into the following proposition for later reference 

Proposition 2.7 The matrix kernel M n (z,w) 

PnW*,«0 := r£qf\w)pM{z) + ¥j ^ 1 ) (2-29) 

is given in terms of the solution of the Riemann-Hilbert problem (2-24) - (2-22) as 

M. n (z,w) :— (2-30) 

z + w 

The relevance of the matrix kernel H n (z, w) will become clear when we will discuss the spectral statistics 

of the matrix model. 

3 Matrix Models 

Consider the vector space TCn of Hermitean matrices of size N endowed with the U(N )-invariant Lebesgue 
measure 

dAf := J] d5R(M iJ )d3(M u ) JJ dM« . (3-1) 

i<j 3 

Inside Hn we consider the convex cone of positive definite matrices 

H+ := {M = M\M > 0} (3-2) 

with the induced measure. Let a, /3 be two positive densities of finite mass supported on the positive real 
axis: a(M) will simply mean the product measure on the eigenvalues of M. 
Define now the measure on Ti~^ x H^ = {(Mi,M 2 )} 

a rut /^f^ am am a ( M i)P(M 2 ) 

MM, ,M 2 ):= dM ld M 2 det(Mi + M2)jv 0-3) 

Definition 3.1 For the positive finite-mass measure (3-3) we define the partition function as the 

integral 

Z N := ( d/i(M u M 2 ) (3-4) 



H+XH+ 



10 



The resulting random matrix model falls into the general class of two ~matrix~mo dels although the 
"coupling" term is not the most common one (which is e Tr ( MlM2 )). However these and much more general 
couplings have been considered relatively recently in [25]. 

As customary, we re-express the Lebesgue measures dMi in terms of the normalized Haar measure 
dUi of the unitary group U(N) and the Lebesgue measure on the cone M. 1 ? 

dMi = G N dUj[dXjA 2 (X) , dM 2 = G N dV J[ dy 3 A 2 (Y) , (3-5) 

where A(X) := det[a^~ ] is the Vandermonde determinant associated with the n-tuple X of ordered 
eigenvalues x\ < . . . < x n of Mi and A(Y) is defined in the same way relative to the matrix M^. The 
constant Gm is not of much relevance. It depends only on N but not on the densities. The resulting 
measure involves two copies of the unitary group and one of the integrals can be performed leading (up 
to a multiplicative constant depending on the normalizations of the Haar measures and the size N) to 
the measure below on the spectra of Mi , M% 



d^(X,Y) = G 2 N A 2 (X)A 2 (Y)( f 

\JUi 



<lL a(X)dXp(Y)dY (3-6) 



l u{N) det(X + UYW) N 

N N N N 

a(X) := J] a( Xj ) , dX := JJ dx 3 /3(Y) := J[ (3( Vj ) , dY := JJ dy 3 . (3-7) 

j=i i= l i= l i =1 

At this point we need to compute the integral over U(N); we can use the result on pages 23-24 of [25] 
for the special example (A-28) 

det \(l-za i b i ) N - r - 1 ] 1 .. .,.. 
det(l -zAUBU^dU = C N , r ^ M lL, m ' "' (3-8) 

IU(N) &{ A )&{B) 

Uk=i k - 

CN ' r '' = z^-^Uk=i(r-N+l) k 

where 

k 

(r-JV + l)*:=JJ(r-JV + j) (3-10) 

is the Pochammer symbol. 

Setting z = -l, A = X' 1 , B = Y in (3-8) we obtain 



dU det 



(1 + ^)iV-r-l 
_ p L ' J l<i,j<N (''?-11 s | 

u{N) det{l + X- l UYW) r N ' r A(X-!)A(F) { ' 



Multiplying both sides by det(X) r we then obtain 



L 



dU = {-IFWC*,** [iXi V/^7l ]l ^^ (3-12) 



u{N) det(X + UYWY y ' x ' A(X)A(Y) 

11 



The case of main relevance to us is r = N, which yields (Cn,n = (—1) N ( N 1 ^ 2 ) 



/ 



dU 



det 



Hlj 



(3-13) 



l u{N) det(X + UYW) N A(X)A(Y) 

This shows that the measure d/x(Mi, M 2 ) can be reduced to a measure on the spectra of the two matrices 
(we use the same symbol for the measure on the eigenvalues) 



dn(X,Y) = G%A 2 (X)A 2 (Y) de ^[ X ^ ] a(X)dXf3(Y)dY 



A(X)A(Y) 



K(x,y) 



x + y 



(3-14) 



In fact we could have used the general formula (3-8) for any r := N + h; note, however, that for h 
integer and less than —1 Harnad-Orlov's formula in the form presented above cannot be used, since the 
determinant in the numerator vanishes (for N > —h) and so do some denominators in the definition of 
the constants Cjv, r - In other words, one should take an appropriate limit and use de l'Hopital's rule. For 
any h > or any h $■ — N, tracing the steps above one could obtain general models where the reduced 
measures on the spectra have a form 

det (K h (xi,yj)) 



d» h (X,Y) = G' N C N , r A\X)A\Y)- 



corresponding to the unreduced measure 



A(X)A(Y) 

K h (x,y) 



a{X)dX(3{Y)dY 
1 



(x + y) 1+h 



dnh(M 1 ,M 2 ) = dMidM 2 - 



a(M 1 )(3(M 2 ) 



(3-15) 



(3-16) 



' det (Mi + M 2 ) N + h ■ 

It will be important for us in what follows that for any value of h ^ — N the kernel K^ is totally 
positive or at least sign regular on M + x R + . (Sign- regularity means that determinants in Def. 2.1 are 
nonzero and their sign depends only on their size.) 



Lemma 3.1 Consider the kernel Kv, 



(x+y) 1 + h ' 

1. for 1 + h > 0, Kh is totally positive, 

2. for 1 + h £ R_ \ — N, Kh is sign-regular. 



restricted to (x, y) £ M+ x 



Then, 



0, 



. , n, s 



Proof. To prove this assertion define for < yo < • • ■ < y„ functions ki{x) = , ,. , 

1 + h . Denote for j = 1, 2, ..., Cj(s) = s(s + 1) . . . (s + j — 1) and Co(s) = 1. Then the Wronskian 

W(ko, . . . , k n )(x) of ko, ■ • • , k n is equal to 



det 



(-iy Cj (s) 



1 



(x + yiY+i 



= (-1) 



1,3=0 



TT^Vdet 



1 
(x + y t y 



(3-17) 



i,j=0 
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Thus, W(ko, . . . , k n )(x) is a nonzero multiple of (— 1) 2 A(zq, . . . , z n ), where Zi — -7—^ — c and A(z , • • • , z n ) 
is the Vandermonde determinant constructed out of the z^s. Since zq > . . . > z n , we see that 

(-l) M ^ i A( 2o ,...,z„) = A(z„,...,z ) >0 (3-18) 

Thus W(ko, . . . , k n )(x) is positive for any n if S > since all the Cj(s) are positive numbers; if s < 
then - denoting by [s] the greatest integer less than s (hence negative) - we see that the sign of W is 

sign{W) = (-\)-^ n - ljM i M1 (3-19) 

for any n > and x > 0. Together these observations imply (Thm 2.3, ch. 2 of [31]), that Kh is totally 
positive for s > on R + x K + or at least sign regular for s <G R_ \ — N. Q.E.D. 

In App. A we will show that for h integer or half-integer, the model (3-16) is the reduction of a 
3-matrix model, where a "ghost" gaussian matrix A has been integrated out. Depending on the range 
h < N or h > N the matrix A consists of ordinary variables (bosons) or Grassmann variables. 

Remark 3.1 The fact that the measure dfih depends on N (in the determinant in the denominator) is 
a feature of the model rather than a problem; indeed, in studying the large N limit it is natural to make 
the strength of the interaction increase at the same rate as the size of the matrices. For example, in the 
"standard" two-matrix model one considers the interaction e cNTr ( M i M 2) 

3.1 Correlation functions: Christoffel— Darboux kernels 

In this section we will compute the correlation functions of the model. More precisely, 

Definition 3.2 The correlation functions of the model are defined as 
TZ {r ' k) (x 1 ,...x r ;y 1 ,...,y k ) := 



(N - r)\(N - k)\Z N J 11 y J> 11 A(X)A(Y) 

where Z N := ^ J A(X)A(Y)det[K(x i ,y j )] i , j < N a(X)l3{Y)dXdY. 



(3-20) 



These functions allow one to compute the probability of having r eigenvalues of the first matrix and k 
eigenvalues of the second matrix in measurable sets of the real axis. 

The computations of these correlation functions in term of biorthogonal functions reported below 
follows the general approach in [19] and [23]. 

Using the well-known formula for the Cauchy determinant 

1 

. X i+Vj. 



det 
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we obtain the measure 

A(A) A(Y) ( X )p( Y )dXdY (3-22) 

Uiji^ + Vi) 

We will be using the correlation functions only to compute expectations of spectral functions, namely 
functions of A, Y which are separately symmetric in the permutations of the x/s or y/s. 

Lemma 3.2 Suppose F(X) is a symmetric function under the action of the symmetric group Sn- Then 

i /■ r N i 

— / F(X)A(X)det[K(x j ,y j ))a(X)dX = / F(X)A(X) J[ ——a{X)dX 

■ •* J 7 = 1 ^ ^ 

Proof. Under the assumption F(X) = F(X a ) for any a £ Sn we have 

N 

/ ffflAffl II —-a(X)dX = 

J 7=1 j ' y j 

- M E J nX)A(X a ) n — - — a(X)dX = 

^ [f(X)A(X) ]T e(a)f[K(x aU) , yj )\a(X)dX^ 



TV! 

^y / R .V)A( A') dot [A'(,-,. /,;)],,(. YM.V 



Q.E.D. 



Corollary 3.1 Lei -F(A, Y") &e symmetric with respect to either set of variables X or Y. Then 

N 

^ J F(X, Y)A(X)A(Y) det[K(x j ,y j )]a(X)0(Y)dXdY = J F(X, Y)A(X)A(Y) f[ _l_a(X)/J(Y> 

(3-23) 



Since we are interested in the unordered spectrum of the matrices M\, M2, in view of the above Cor. 3.1 
we will focus henceforth on the following unnormalized measure 

N 

/ = i x j + Vj 

From the properties of the Vandermonde determinant we can write 



dV{X,Y) ■- A(A)A(F) J| a ^ Xj)l3{ - V]) dXdY = A{X)A{Y)-^^-dXdY (3-24) 



A(A)A(F) = det \pi( Xj )} det [§•(%)] (3-25) 

where pj, qi are any monic polynomials of exact degree i in the respective variables (i = 0, . . . , iV — 1). 
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It is therefore natural [19] to choose the sets of monic polynomials {pj(x),qj(y)}j^ to be biorthog- 
onal with respect to the bi-measure " ^1T clxdy. This means that 

/ p k (x)Uy f {x) ! 3{y) dxdy = c k 2 5 ke . (3-26) 

s + JR+ x + V 

The constant c k was defined in (2-4, with K(x,y) = l/(x + y)). The normalization constant for dv is 
thus 



-■N 



JjMX,Y) = M ( n cU = Nldet [l^] < J . fc < JV _ 1 (3-27) 



If we introduce the orthonormal biorthogonal polynomials 



p n := —p n , q n := — g„ (3-28) 



then we can write the normalized measure as 

N 



u(X, Y) = ±- det fa{ Xj )] det My,)] f[ a{Xj) f {yj) d N Xd N Y. (3-29) 

Notice now that the product of determinants in (3-29) is a determinant of the product of the two 
matrices indicated; hence 



det \pi-i(xj)] id < N det fe-i(y»)]ij-<jv = dctK 7V (a; l ,2/ i ) (3-30) 

Kjv(a;,2/) := 53pj(a0%-(y) (3-31) 



The kernel Kjv(a;,y) is a "reproducing" kernel 



I1 A (,-.,/) : // K N (x,z)K N (w,y) a{w)f3< f )dwdZ (3-32) 

2 + 10 



which follows immediately from the biorthogonality In addition we have 

a(x)/3(y)dxdy 



N (x, y) v '™' "- = N . (3-33) 

x + y 



Summarizing, in the new notation we have 

Proposition 3.1 The probability measure on X x Y induced by (3-3) is given by: 

1 



( x]]2 det[K N (x i ,y j )}det[K(x i ,y j )]a(X)P(Y)dXdY (3-34) 

while the correlation functions are: 

TZ^(x 1 ,...x r ;y 1 ,...,y s )= (3-35) 

rn=io(«i)nLij9(«j) r A A 

(N-rV.(N-s)\ det[K N (x l ,y j )]det[K{x I ,y 3 )] [[ a{x 3 )dx 3 [[ (3{y 3 )dy 3 (3-36) 

^ '*'' '' J £=r+l j=s+l 
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Example 3.1 Consider r = l,s = 0. In this case, with the help of Lemma 3.2, we obtain: 

1 /• N N 

^ (1 ' 0) - (N _ 1)[m a ^ / det[K N (x i ,y j )]det[K(x i ,y j )}l[a(x i )dxel[p(y j )dy j = 

1 - N N N 

— — — a(xi) / det\K N (xi,yj)] J[ J[ a{x t )dx t J[ P{y 3 )dy 3 = 

\ !■ J i=1 i "r Vi e=2 3=1 

—— a(»i) ^ Z{v)t{v') I Pa{l)-l{Xl)---P*{N)-l{XN) 



(N 

cr,a'£S N 



N , N N 



K N (xi,yi) 



q*>(i)-i(yi) ■ ■ ■ q<r'(N)-i(yN) Yl — t — II a ( x e.) dx i IJ PivMvj = a ( x i) / - "2'" 1 ' P(.yi)^-yi 

i=l Xi + y * £=2 j=l J Xl + yi 

(3-37) 

3.1.1 Correlation functions in terms of biorthogonal polynomials 

In the paper by Eynard and Merita [19] (generalized in [23]) they never used any specific information 
about the model they were considering (with the Itzykson-Zuber interaction) but only the fact that the 
matrices were coupled in a chain. We recall the relevant result here 4 . Define 

Hao(x,y) :=K N (x,y) 

H 01 (x,x') := f Hm (x,y)^y , H w (y,y') := [ H 00 (x,y')^^ 

J x' + y J x + y 

rr / n ff TT / -. a(z)dzd(w)dw 1 ,„ „„. 

H ll{ y,x) := JJ H w {zM ^l y) ll w) - — y • (3-38) 

Since Hqq is a reproducing kernel 

J H 00 (x, y)H xl {y, x')[3(y)dy = H 01 (x, x') - H 01 (x, x') = (3-39) 

and similarly 

/ H n (x, z)H QQ (z, x')a{z)dz = H w (x, x') - H w (x, x') = . (3-40) 

Integrating these two equations against (x + y')~ 1 a(x)dx (the first) or (x' + y')^ 1 (3(y')dy' (the second) 
we find also 

' H 11 (y,x)H 10 (y,y')l3(y)dy= f [ H 11 (x,y)H 01 (y,y')a(y)dy = (3-41) 



The correlation functions for r eigenvalues x\, . . . , x r of M\ and s eigenvalues t/i, . . . , y s of M% were 
computed in [19] and are given by 

r s 

TZ {r ' s) (x 1 ,...,x r ;y 1 ,...,y s ) = J| a(xj) Y[(3(y k )x 

3 = 1 k=\ 



4 It could be extended to the chain of matrices (see below) but it becomes a bit cumbersome to describe. 

1G 



x dct 



[ H 0l^i^j)]i< itj < r 


[ H w( X i>yi)]l<i<r,l<j< e 


_ [HuiyuXj)]^.^.^ 


[ H w(.V*>Vj)]l<ij<. . 



(3-42) 



where 0<r<N,0<s<N, l<r + s, with the understanding that if cither r or s is then 
the corresponding blocks labeled by r, s respectively, are absent. Thus, for example 7?.( 1,0 )(a;i) = 
a(xi)H ,i(xi,Xi). 

This is a nontrivial result and perhaps one can get a bit of an insight by considering the special case 
r = s = N. We know that in this case 



fc( N ' N) (xi, ...x N ;yi,...,y N ) = 
de%[K N { Xi , Vi )\ deb[K(xi, y )]a{X)(3{Y) 

so we have the following not at all obvious identity: 
Lemma 3.3 



det[K N (xi,yj)} det[K(xi,yj)] = dct 



(3-43) 
(3-44) 



" [Hoi(xi,xj)] itJ < N 


[HoQ(xi,yj)].^ N ' 


_ [fl'ii(w,ajj)] w < w 


[H 10 (y l ,y J )] i ^ N 



(3-45) 



Proof It suffices to observe that the 2N x 2N matrix on the right hand side has the following block 
structure: 



K 



-P (x) J Q 1 (-x) 



P (x) J Q (y) 



(3-46) 



_Pi(-yr Qi(- X ) - K ( x -y) -Pi(-y)Qo(y)_ 

where Po(x) = \pi-i(xj)]i<i,j<N-, Pi(y) = \Pi-i(Vj)]i<i,j<N and, after exchanging p { s with %s, we 
define the remaining symbols (see equations (2-11) and (2-10) for definitions) accordingly. Moreover, 
K(x, y) = [K(xi, yj]i<i,j<N- With this notation in place it is clear that K, admits the following (Bruhat) 
decomposition : 



K = 



Po(x) J 

-Pi(-y) T i 



I 

1 



-K(x,y) 

-Qi(-x) Q (y) 



(3-47) 



Thus det/C = det P (x) T detK(x,y) det Q (y) = det[K N (x i ,y j )]det[K(x i ,y j )]. Q.E.D. 

A comparison of the definitions (3-38) with the entries of the Christoffcl-Darboux identities of Prop. 
2.3 shows that they are intimately related, in fact 

Ott,i<V: 



H^(x,y) = (-r +v Yl qf((-Yx)pt\(-Yy) 



3=0 



x + y 



(3-48) 



and we recognize that (up to some signs) these are precisely the entries of the H„ kernel of Prop. 2.7: 

H^x, y) = {-Y +v W n ^ v ((-rx, {-) v y) (3-49) 

More explicitly we have 
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Proposition 3.2 The kernels of the correlation functions are given in terms of the solution of the RHP 
in Prop. 2.5 as follows 



H 00 (x,y) = M. Nfi0 (x,y) = 

H i(y',y) =-H JVj0 i(-y / ,y) = 
H w (x,x') = -Hj^ioO*:, -x') = 



T-\-y)T(x) 



.3.1 



x + y 
f^(- y )f(- y ') 

y' -y 



:S.2 



Y- x {x')Y(x) 



(3-50) 



2.1 



Hn(y,x) 



lN,n(-y,-x) 



X' — X 

f- 1 (x)f(-y) 
x + y 



The potential relevance of these formulae is that it allows one to compute the large N asymptotic be- 
haviour of the correlation functions in terms of the asymptotic behaviour of only 3 consecutive biorthog- 
onal polynomials and auxiliary functions associated with them that enter in the Riemann-Hilbert for- 
mulation of Props. 2.4, 2.5. 

The steepest-descent analysis of these problems which will appear in a forthcoming paper. 

We expect that, after a complete description of the asymptotics of the BOPs is obtained, these formulae 
can be used to address the issue of universality for this matrix model via the Riemann-Hilbert approach. 

Remark 3.2 Note also that the formula above is not symmetric inasmuch as the BOPs p, q play different 
roles in the Christoffell-Darboux theorem; we could however rewrite the theorem as 

JV-l 

(3-51) 



^2pj(%)qj(y) 

3=0 



PI N (x)B N (-x)ci 0!N (y) 



x + y 



where B(z) was defined in (2-9). The auxiliary vector q 0l iv enters in a similar Riemann-Hilbert formu- 
lation with the roles of the densities a and [3 interchanged. 

3.2 A multi matrix model 

It is possible to extend the model (3-3) to a chain of matrices with the nearest neighbor interaction 

' (3-52) 



N+hi 



det(Mi + Mi+i) 

The "strength" of the interaction (i.e. hj's) may depend on a site along the chain. Specifically, consider 

the space (Ti^r) an< ^ a collection of R positive densities ct\,.. . ,a_R on M + . Define the finite mass 

measure 

1 



n-i 



R 



d»(M 1 ,...,M R )= Yl 



det(M e + M e+1 ) N + h ' 



^Y[a e (M e )dM e 



(3-53) 
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Using the Harnad-Orlov formula one we obtain, up to normalization constant, the following measure 
which we denote by the same symbol 

JL ^}det\(x i i+Xi +li )- 1 - he ]^. .... 

MXi, ...,X R ) = J] a,(X,)A(X,) 2 dX, p \ { x e )A(X e+1 ) (3 " 54) 

Following the same steps that led to the expression (3-24) as positive density for (S'jv) 2 -invariant 

observables, we obtain the following reduced measure on the spectra (up to a normalization depending 

only on /i/s but not on the measures) 

fl-i R 

dv{X 1 ,...,X R ) = A{X 1 )A{X R ) I]det(X £ + X £+1 )- 1 -^ Y[a e (X e )dX t (3-55) 

£=1 1=1 

In the case where all interactions are the same hi = we have 

R-l R 

dV{X t , ...,X R ) = A{X 1 )A{X R ) J] Y[a e (X e )dX e (3-56) 

1=1 \ £ l+l ' 1=1 

which can be seen as a generalization of (3-24) . This model too can be treated with the aid of biorthogonal 
polynomials that will now satisfy a R + 3 recurrence relation; note that the length of the recurrence 
relations does not depend on the "potentials", or densities, Oj. This is in sharp contrast with the usual 
multi-matrix model with interaction Q cTc i M i M i+^) [19, 6]. To characterize these biorthogonal polynomials, 
an (R + 2) x (R + 2) Riemann-Hilbert problem can be set up and the strong asymptotics can be dealt 
with, but the complexity of this problem definitely warrants a separate paper. 

4 Diagrammatic expansion 

In parallel with the -^ expansion for the Hermitean matrix model and the IZHC two-matrix model, we 
would like to sketch the similar formal expansion of the model in terms of colored ribbon graphs. The 
weights a and (3 entering the definition (3-3) are assumed to be of the form a{M\) — Q- NTr ( u ( M ^)> 5 f3{M2) = 
e -jVTr(y(Mi)) _ y^ e p er f orm a shift and a rescaling of the matrices so that det(Mi + M%) *-* det(l — C,M\ — 
77M2). Of course the values of C, and 77 are suitably restricted to a neighborhood of the origin: however 
this restriction is irrelevant since the manipulations below are in the sense of formal power series. 

The procedure amounts to a perturbative Taylor-expansion around a Gaussian integral. In other 
words we will be considering a partition function in the form 

Z N := /7dAfidA/ 2 e _ ^ Tr ( Ml2+M22 ) +WE ^ lTr(CMl+ ^ M2) ^ WTrap(Ml)_WTry! ' (M2) . (4-1) 



where U P ,V P are the perturbations of the Gaussian (including a quadratic term as well) which, for 
convenience, we parametrize as the following formal series 



3C 



U p (x) := - £ ^4^ , V p (y) :=-J2 *V^ ^ 
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(4-3) 



where the average is taken w.r.t. the underlying {uncoupled!) Gaussian measure 

N 



-y— dMidM 2 exp 



-yTr(M? + Aff) 



Z ( ° ] := j dMidM 2 exp 



-'—Tr(M? + M2) 



(4-4) 



Note that, due to the shifts in Uj,Vj, the £ — 1 term in the third sum above cancels exactly against the 

nn \\r 4- in o farm 



shifts in the first two sums; similarly, only the term ^Tr(MiM 2 ) remains in the quadratic part: 



Z N = {Jexp NTi (uiMx + y Ah 2 + v x M 2 + y M 2 2 + C,r)M x M 2 + •■•)) (4-5) 

Using Wick's theorem for the evaluation of Gaussian integrals and the frequently used combinatorial 
interpretation (see [15] for an excellent introduction) one sees that the partition function is the sum over 
all possible bi-colored ribbon graphs respecting rules which we now specify 

• There are two colors for vertices/edges (say, red/blue); vertices are distinguished in monochromatic 
and bichromatic. 

• Each monochromatic vertex of valency j enters with a weight Nuj/j (blue) or Nvj/j (red). 

• Each edge (red or blue) enters with a weight -^ . 



k„€ — k 



• Each bichromatic vertex of valency £ enters with a weight ^ — where l<fc<^ — lis the number 
of blue half-edges and t — k the number of red ones and appears in with a multiplicity I k ) corre- 
sponding to all possible arrangements of k blue legs amongst £, up to cyclic reordering. In particular 
there is only one bichromatic bivalent vertex (up to automorphism) which enters with weight (f]/N. 
In general, since each such vertex corresponds to a trace of the form Tr(Mi ai M 2 1 Mi°' 2 . . .) and 
in view of the cyclicity of the trace, there are precisely £ equivalent vertices obtained by cyclically 
permuting the matrices in the sum, which corresponds diagrammatically to a rotation of the colors 
of the legs of the vertex. Hence there are in fact \ ( k ) = (i-kYk\ mec l u i va l en t bicolored ^-valent 
vertices in each diagram contributing with a weight ( k rj i ~ k /N. 



• 



Each connected Feynman diagram T constributing to the perturbative sum has a power N 
N 2 ~ 2g where g = g(T) is the genus of the surface over which the graph can be drawn. 
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Summing up over all possible labeling of the bicolored fat-graphs 
r leaves a factor |ylu£(T)| in the partition function (see [15]) and the 
result is hence 



InZ 



A' 



jy2-2 gr 

^ \Aut(T)\ 
Connected J- 1 

bichromatic 

fatgraphs F 



iks nnc 



(i-Qri^kri 



(4-6) 



=2fe=l 



where rij = Tlj(T) is the number of blue vertices with valency j, rrij = 
mj(T) the number of j-valent red vertices and r k = r e k (T) is the number 
of ^-valent bichromatic vertices with k red (and hence £ — k blue) legs. 



5 Large N behaviour 

Consider two densities a(x),/3(y) of the form 

a(x) = a h (x) = e-^W , p( y ) = p K ( y ) = e -^(y) ? 




Figure 1: An example of a con- 
nected diagram contributing to 
rp the partition function. 

h =N> T >°- 
(5-1) 

Here we have introduced a dependence on the small parameter h on the 

measures (but we will not emphasize this dependence in the notation). From the experience amassed in 

the literature on the ordinary orthogonal polynomials we start by considering the heuristic "saddle-point" 

for the partition function Zn, namely the total mass of the reduced measure (3-24). The heuristics calls 

for a scaling approach where we send the size of the matrices N to infinity and the scaling parameter H 

to zero as 0(1/N), namely that we vary the densities a, j3 by raising them to the power i. We thus have 

1 



Z N = / A(X)A(Y)det 



Vj 



a(X)l3(Y)dXdY 



i,j<N 



Using 3-22 we have 



Zjv oc 



A(X) 2 A(Y) 2 

UiAxi + yj) 



exp 



1 N 

-J2U( Xi ) + V( yj ) 



dXdY 



(5-2) 



(5-3) 



Taking -^ times the logarithm of the integrand we have the expression 

1 N 1 1 1 



TN ^~ w/ ' ' w/ N 2 

3=1 JT^k 

Su{X) + S v (Y) + J(X,Y). 



j^k 



Vk\ 



3,k 



It is convenient -in order to deal with a more standard potential-theoretic problem- to map Y 
and define V*(y) — V(—y) so that we can rewrite the action as 

S{X, Y) = S V {X) + S v * (7) + ^ In |. 



Vk\ 



(5-4) 
-> -Y 

(5-5) 



J-k 
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which describes the energy of a gas of 2N particles of charge +1 (for the Xj's) and —1 (for the y/s) 

separated by an impenetrable, electrically neutral, partition confining the positively charged particles to 

the positive real axis under the potential U and the negatively charged particles to the negative axis 

under the potential V*. 

The usual argument is that the configuration of the minimum contribute the most to the integral and 

-under suitable assumptions of regularity for the potential- one wants to show that the sequence of these 

minima configurations tends in-measure to some probability distributions. 

We define 

1 N 1 N 

P = jy l2 6 *>* > ^ = n X, 6 vi ( 5 - 6 ) 

fc=i j=i 

and rewrite the action as 

S[p,(j\:—— J U(x)p(x)dx — J J p(x)p(x')ln\x — x'\dxdx'- 

1 

+ f 



R+ 

V*(y)fi(y)dy- I I n{y)n{y')hi\y - y'\dydy'- 



P{x)p{y) In \x - y\dxdy (5-7) 

5.1 Continuum version: cubic spectral curve and solution of the potential 
problem 

We immediately rephrase the above minimization problem in a continuum version. In order to emphasize 
the symmetry of the problem it is convenient to denote U(x) by V\(x) and V*(y) — V(—y) by V 2 (y) and 
denote the corresponding equilibrium densities by pi(x), p2(y)- 

The main point of this section is Thm. 5.1 which states that the resolvents (Markov functions) of the 
equilibrium distributions are solutions of a cubic equation that defines a trigonal curve 5 ; this result is the 
analogue of the better-known result for the equilibrium measure appearing in the one-matrix model ([12] 
and references therein) . The derivation that we present here is of a formal heuristic nature inasmuch as 
we discount several important issues about the regularity of the equilibrium measures. However this sort 
of manipulations is quite common and they can be obtained also from the loop equations as done in [18] 
for the 0(n) model. 

Rewriting the functional in the new notation for the potentials we find 

S{pi,p 2 ] := — / Vi(x)pi(x)dx - / / p 1 (x)pi(x')ln\x - x'\dxdx'+ 



Tj R+ 
■fy / V 2 (y)p 2 (y)dy - J J p 2 (y)p2(y')ln\y-y'\dydy'- 



3 Namely a curve of the general form w s + Aw + Bw + C = 0, with A, B, C smooth functions of the spectral parameter 
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Pi(x)p 2 (y)ln\x-y\dxdy (5-8) 

We will be studying the measures p\ , p 2 that minimize the above functional; here we will assume their 
regularity (which shall be proved in a separate paper) . Also, the potentials V\ , V 2 will be assumed real 
analytic. 

More precisely, we make the following assumption on the nature of the potentials 

Assumption 5.1 The two potentials Vi(x),V 2 {y) are restrictions to the positive/negative axis of two 
real-analytic functions such that 

• lim^+oo ^2 = +oo, and similarly lim^^ j|j^ = +oo; 

• the derivatives VUx) have at most finitely many poles in a strip of finite width around R; 

• finally, V\(x) > Ci|ln(x)| for some constant C\ > and x > and V 2 (y) > C2 j In | J/ 1 | for y < 
and some constant C2 > 0. 

It will be shown in a separate paper that this assumption is sufficient to guarantee that the equilibrium 
densities exist, have compact support, are regular and their supports do not include the origin. In order to 
enforce the normalization of pi, p 2 -as customary- we will introduce in the action two suitable Lagrange 
multipliers 

S ■= S[7i,72, Pi, P2] := S[p 1 ,p 2 ] +7+ f / Pidx- lj +7_ ( / p 2 dx- l) (5-9) 

The vanishing of the first variation of the functional 5* yields the equations 

2 / pi In \z — x\dx + I p2 In \z — x\dx = 7+ , z € Supp(pi) 

JR4. J~R- 



T 



+ 
2 / P2 In \z — x\dx + I p\ In \z — x\dx — 7_ , zG Supp(p2) (5-10) 



T 
Differentiating w.r.t. z yields 

}-V{{z)-2P.V. [ pi-*— dx+ [ p 2 - i —dx = 0, ze Supp( Pl ) 
1 J z — x J z — X 

%V 2 \z)-2P.V. f p 2 ^— dx+ I Pi-*— ds-0 , ze Supp{p 2 ) (5-11) 

1 J z — X J z — X 

where P.V. J indicates the Cauchy principal value. These equations are best written in terms of the 
resolvents (or Weyl functions) 

Wi {z) := / M^, zeC \ Supp( Pl ) . (5-12) 
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Indeed, with the help of the Sokhotskyi-Plemelj formula, the equations above take a simpler form 

W ll+ (z) + Wi,_(«) = ±V{[z) + W 2 {z) , z G Supp(pi) ,* ^ 

W 2 ,+ {z) + Wi,.{z) = iV 2 '( X ) + Wi(«) , * G Supp(p 2 ) l °" 1,3J 

Wi,+ - W it - = -2^ (5-14) 

Note that in the RHS of the system of equations, the resolvent of the other measure is evaluated at a 
regular point due to the fact that the two supports are disjoint. 
We define, partly motivated by [21], the shifted resolvents 

2V1 + Vi V7 + 2K' 

Yi := -W 1 + ^ 2 , F 2 ~ VK 2 - l J T 2 . (5-15) 

Observe that, in view of Assumption (5.1) , {Yi, Y2} have the same analytic structure as {W\, W2}, while 
the equations describing the jumps simplify to: 



ii,+ + *i,- = -Y 2 (z) , z e Supp{ Pl ) 
Y 2 ,+ + Y 2 ,_ = -Yi(z) , z € Supp{p 2 ) 



(5-16) 



*i,+ - *i,- - 2i^ Pl , y 2 ,+ - y 2 ,_ = -2iTr P2 (5-17) 

Using these equations one obtains by direct inspection that 

R(x):=Y 1 2 + Y 2 2 + Y 1 Y 2 (5-18) 

has no jumps on either supports. Multiplying eq. (5-18) on both sides by Y\ — Y 2 one obtains 

R(z)(Y 1 -Y 2 ) = Y 1 3 -Y 2 3 (5-19) 

which can be rewritten as 

Yx 3 - R(z)Y 1 = Y 2 3 - R{z)Y 2 := D(z) . (5-20) 

The first expression may a priori have at most jumps across the support of pi, while the second may 
have jumps only across the support of p 2 : since the two supports are disjoint we conclude that D(x) is a 
regular function on the real axis. 

Remark 5.1 If we introduce the function Yq := —Y\~Y 2 the jump relations (5-16, 5-17) can be rewritten 
(after a few straightforward manipulations) as 

Yq.±(z) = Yl t (z) z e Supp(p!) , , 

Y ,±(z) = Y 2 ^{z) z e Supp( P2 ) . [ °- ZL > 

This implies that we can think of the Yj(z) 's as the three branches (three sheets,) of a Riemann surface 
realized as a triple cover of the x-plane branched at the endpoints of the spectral bands. The Riemann 
surface is depicted in Fig. 2. Thm. 5.1 will be proving more formally this statement, realizing Yq,Yi,Y2 
as branches of a cubic equation (which the reader can already guess from (5-20)). 
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Figure 2: A pictorial example of the three sheets of the Riemann surface of Y(x). In blue the support of 
pi (right) and in green that of p 2 (left) . 



Theorem 5.1 Under Assumption (5.1), the two shifted resolvents satisfy the same cubic equation in 
the form 

E{y,z):=y z -R{z)y-D{z) = Q . 



The coefficients R(x), D(x) are related to the equilibrium measures as follows: 
D(x) = UqUiUi + U 2 Ri + U 1 R 2 + iViaK - ^V 2 i^ 



(5-22) 



(5-23) 

(5-24) 
(5-25) 



where the functions Uo,U\,U2 and the operators Vi, V2, Vi2t V21 have been defined as 



Ui:= ?K±K, U2 



3T 



V,/(z) := 
V 12 /(z) := 

V 21 /(z) := 



/(*) - /(*) 



2V 2 ' + v/ 

3T 
dft(aj) , 



, ^o 



3T 



z — x 
f(z) - f(x) dpi(x)dp 2 (y) 

z — x x — y 

f(z) - f(x) dp 2 (y)dpi{x) 

z-y y-x 



(5-26) 

(5-27) 
(5-28) 



Furthermore the support of pi,p 2 must consist of finite union of finite disjoint intervals. 
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Proof. The statement about the support of the measures follows from the argument used in [12] 
that that we sketch here. By results of Chapter 13 of [36]), Assumption 5.1 implies that the supports of 
the equilibrium measures are compact. Moreover, it follows that the branchpoints (and branchcuts) of 
the functions Yi,2,o coincide with the supports of pi,p%. Since they solve a cubic algebraic equation, the 
branchpoints are determined as zeroes of the discriminant A := 4i? 3 — 27 D 2 of (5-22). 

As suggested by the previous discussion, even without an explicit expression for R, D, it follows from 
Morera's theorem that R, D are also real-analytic and so is A. Thus A cannot have infinitely many zeroes 
in a compact domain and so the number of endpoints of the supports is a-priori finite. 

The only part that is left to be proven are the formulae for R, D. 

Computing the jump of Wj from eqs. (5-14) one obtains 

(Wi)+ - (Wi)l = -2itt Pi (hr{ + W 2 ) (5-29) 



T 
and a similar equation for W 2 . This implies that 

Wi2 rp_Momii±^m dx= 

I z — - X 

1 v , ( s w ( , f Pi{x){V{{z)~V{{x)) A ff Pl (x)p 2 (y)dxdy 

= -v l{z )w l{z) -J ^-^ d * + JJ {z - x){x - y) ( 5 " 3 °) 

=:fli(z) 

Note that R\(z) is regular on the support of p\ since V{ is. Hence we have 



Wl \ z) = l -v[{ Z )w l{z) - Rl{Z ) + // f;! p ^ d _ x y" ( 5 - 31 ) 



2i..\ _ l irt/~\txr /~\ d /„\ , ff Pi(x)p2(y)dxdy 



wAz) = f vi( z) w 2 ( z ) - M z) + JJ jt^W ■ (5 " 32) 



:W 21 



Adding the two together and using the identity 

1 



(5-33) 



(z-x)(x-y) (z-y)(y-x) (z-x)(z-y) 

we obtain 

Wi 2 + W 2 2 = yV{Wx + yV£W 2 -Ri-R 2 + W X W 2 , (5-34) 

which is precisely eq. (5-18) when rewritten with W\ 1 W 2 replaced by their expressions in terms of the 
Y\,Y 2 variables resulting from (5-15). 

The second formula for D(x) can be obtained by noticing that D(x) = Y^Y\Y 2 . Explicitly, this reads 
(after rearranging the terms and using the definitions for the shifted resolvents (5-15) and the Uj'a (5-26) 

/ V' \ / V 

^^Jir 1 11/ TI/2 tt I ti/2 1 iir \ TT 1*^2 2 1 



D(x) = U UiU 2 - 2U W 1 W 2 - W(W 2 + WiW$ - U 2 I W( - -^W x \-U x \W*- -±W% ) (5-35) 
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Substituting in the last two terms the expressions (5-31, 5-32) we obtain 

D(x) = U0U1U2 - 2U WiW 2 - WfW 2 + W{Wl + U 2 Ri + U-lR.2 - U 2 W 12 - U x W 2 i (5-36) 

Using W{W 2 = W 12 + W 21 (5-31, 5-32) we obtain 

D{x) = U U X U 2 - WfW 2 + W x Wl + U 2 R X + U X R 2 + ^W 12 - ^W 21 (5-37) 



Now we can rewrite 



v;w 12 = v 12( y/) + // y((fp^(y)^y (5 _ 38) 

JJ {z-x){x-y) 
V 2 W 21 = V 21 (V 2 ) + // _ (5-39) 



so that we have obtained 



D{x) = UoUiUz + U 2 R 1 + U X R 2 + ^V 12 V{ - ^,V 21 V 2 ' + m (5-40) 

K:= l -[[(YM + YM) eMpM^l _ W?W2 + WlWi (5 _ 41) 

TJJ \z-x z-yj (x-y) 

We want to show that 9^ = 0: indeed it is clear that the only part of D(x) which may have jump- 
discontinuities is £H, but we know that D{x) has no such discontinuities. 

Hence $K has no discontinuities; on the other hand it is clear from its definition that d\ cannot have 
any other singularities and hence it is an entire function. Inspection shows that 9t(z) — > as z — > 00 and 
hence 91 must be identically vanishing by Liouville's theorem. This concludes the proof of the theorem. 
Q.E.D. 

The presence of a "spectral curve" will be one of the crucial ingredients for the large-degree asymptotic 
analysis of the biorthogonal polynomials using the Riemann-Hilbert formulation given in [8] and the 
Dcift Zhou nonlinear steepest descent method. Indeed we will show in a separate publication that the 
OPs are modeled by (spinorial) Baker- Akhiezer vectors (similarly to [9]) that naturally live on the three- 
sheeted covering specified by (5-22) . 

A A rectangular mixed 3— matrix model with ghost fields 

The model (3-16) (for any value of h) can be obtained from the following matrix models 

A.l Integer h < N 

Consider the standard Lebesgue measure on the space Mat((N — h) X N,C) of complex (N — h) X N 
matrices, viewed as a linear space: if A £ Mat((N — h) x N, C) we will use a shorthand notation 

N-h N 

dAdA^ := Y[ Yld^A^d^A.j (A-l) 

»=1 3=1 
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for the volume element. Let, as above, M\,M% £ Ti^ and consider the following normalizable measure 
on the space H% x H% x Mat(N - h, N, C) 

di2c(M 1 ,M 2 ,A) ■- dM 1 dM 2 dAdA t a(M 1 )/3(Af 2 )e- TrA(Ml+M2)At (A-2) 

Since the measure is Gaussian in A one immediately sees that (up to inessential proportionality constants) 

d£ c (Afi , M 2 , A) ex d(i h (M i , M 2 ) 



/ 

■7 i\ 



Mat(N-h,N. 

= dM ld M 2 a(Ml) ^ (M2 



det(Mi + M 2 ) JV - h ' 
/ie{0,l,...,JV-l} (A-3) 



A. 2 Half-integer h < N 



One could also consider a similar model where Mat(N ~ h, N, C) is replaced by Mat(2N — 2h, N, K) 
where h £ \Z and h < N obtaining then 

IN-lh N 

dA-.= n ri dA ^ 

d/I R (M 1 , M 2 , A) := dM 1 dM 2 dAa(M 1 )/3(M 2 )e- Tr ' 4(Ml+M2) ' 4t (A-4) 

/ d-p n (Mi,M 2 , A) ex d/*(Mi, M 2 ) = 

JMat(N-h,N,S) 

= dAf ld Af 2 -4P4^) , (A-5) 

det(Mi + Afa)^' 1 

Of course, if /i is actually an integer then this case reduces to the previous one. 

A. 3 Integer and half integer h > N 

For the sake of completeness we note that when h > N the determinant in the reduced measure (3-16) 
is a positive power in the numerator; this can be obtained from a Gaussian integral over Grassmann 
variables which can be obtained from "complex" anticommuting variables (for h integer) or "real" (for 
h £ |Z) much along the line of the previous (commuting variable) case. 

Clearly, the plethora of models is endless and each value of h can be studied with the aid of a different 
set of biorthogonal polynomial. 

We singled out the case h = 0, since it corresponds to the Cauchy BOPs that have many features in 
common with the much better-known orthogonal polynomials and yet seem to have a very rich but still 
tractable asymptotic theory. 
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B Relation to the 0(1)— model 

The 0(n) matrix model [18] is a multimatrix model for a (positive) Hermitean matrix M and n Hermitean 
matrices Aj with distribution 



AM Y[ dAj exp ( -JVTr 



3=1 



(B-l) 



With some manipulations, the integration over the Gaussian variables Aj can be performed and the result 
written in terms of the eigenvalues Zj of the matrix M, yielding a (unnormalized) measure over the space 
of eigenvalues given by 

A(Z) 2 £ dz lG - NV(z ^ A(Z) 2 a(Z)dZ ,. , . m . 

We wish to show that we can specialize our Cauchy two-matrix-model so that, in principle, it reduces to 
a 0(1) model. More precisely, we will show that the partition function of the Cauchy two-matrix model 
is the square of the partition function of the 0(1) model for a particular choice of measures (see Prop. 
B.l). 

We recall that from a diagrammatic standpoint the large N expansion of the 0(n) model describes 
formally a gas of self-avoiding loops of n different colors (hence in our case of a single color) on random 
surfaces. 

Proposition B.l The partition function of the Cauchy 2-matrix-model with N = 2k and (3{x) — xa(x) 
is the square of the partition function of the 0(1) -model of the same size. More precisely 

( 1 f A(Z)ig(Z)dz \ 2 W J Y\f=i yj A(X)A(Y)a(X)a(Y)dXdY 
\2 k klNlJ K * Ili<j(«i + «*)"/ Wr-A- " dct(X + Y) [ ~> 

Remark B.l We conjecture an analogous statement to hold for odd N. 

Proof. The matrix of moments for this choice of measures reads Iij = jf x t y : ' +l a(x)a(y) £ v and the 
partition function for the corresponding size-TV Cauchy matrix model is by equation (3-27) (up to some 
combinatorial coefficients) the principal minor of this matrix det(/)jv (we use the subscript N to denote 
the principal minor of size N) . Such minor coincides with the RHS of (B-3) 
We introduce the new the skew symmetric matrix 

i rr , ,, . d 2 « 



M ij =-JJxY(y-x) J ^ y (B-4) 

which is clearly skew-symmetric. Denoting by a — [ J x J a(x)dx] t j =0 x the infinite vector of moments of 
the measure a(x)dx, we see that 

I = M+ l -aa l (B-5) 
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and hence 



det(Jjv) = det(Mjv) + -a*Mjv 



(B-6) 



where the tilde denotes the classical adjoint matrix of the principal minor of size N x N (the transposed 
of the matrix of cofactors) . 

Since we are dealing with the case N = 2k even, the adjoint of a skew-symmetric matrix is skew- 
symmetric as well and hence the second term has to vanish. Thus we have 



det(/) JV = detMjy = (Pf(M N )) 2 
We see that up to a multiplicative factor -r^ 



(B-7) 



k\2 k Pi (M N ) = J2 eWjIV^ 
aeSjv 3=1 



a- v 



n 

3=1 



Xj + yj 



1 Vi 

xi v\ 



<-' v N , 



1 Vk 
x k V 2 k 



<-' V? 



We denote by Z the vector of size N with components 

Z = (Z1,Z2,...,Z N ) = (Xl,yi,X2,V2,- ■■ ,Xk,Vk) ■ 

With this notation in place we have 



(B- 



fc!2 fc Pf M N = 



*(z)dzA(z) n - j ^- 

f = \ xj + y 3 



I n — !— n%- n fo+w) n fe+^x^+ww^A^dz (b-9) 

/DiV „ 2l + 2/ 



j<i<N 3 j=l j^£<k j<£<k 

v '" v ' 

T(Z) =:R{Z) 

The rational function T(Z) is invariant under permutations of the variables, whereas R(Z) is not; however 
R(Z) has the same total degree as A(Z), in particular its degree in j/i is N — 1. 
Symmetrizing under the integral sign the variables gives: 



fc!2 fe PfM w = ^ J2 J T(Z (7 )A(Z a )R(Z a )a(Z)dZ = ^ J T(Z)A(Z) £ e(a)R(Z a )a(Z)dZ. 

(B-10) 
Now note that necessarily 

J2 e(a)R(Z a ) = A(Z). (B-ll) 



<j£Sn 
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Indeed R(Z) containts the monomial y Y 1 (with the coefficient 1), and has the same total degree in all 
variables; since the antisymmetrization must be divisible by A(Z) the assertion follows. We thus have 

1 f A 2 (7) 
Pf Mn = »,,,,„ / ^ -r 2 MZ)AZ. (B-12) 

We see then tha up to a proportionality constant Pf Mjq is the total integral of the measure (B-2) 
and thus the proposition is proved. Q.E.D. 

The relationship between the two model does not seem to go much further in the sense that there is 
no direct and simple relationship between the correlation functions of the two models. It seems, however, 
that some connection should be present and is worth exploring. We leave it as an open problem to 
establish a connection between these two models on the level of the correlation functions. 
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